Understanding differential equations through diffusion 
point of view: non-symmetric discrete equations 



o 

u 

Or 

< 



< 
o 



> 

(N 
O 



X 



Dohy Hong 

Alcatel-Lucent Bell Labs 
Route de Villejust 
91620 Nozay, France 

dohy. hong® alcatel-lucent, com 



ABSTRACT 



In this paper, we propose a new adaptation of the D-iteration 
algorithm to numerically solve the differential equations. 
This problem can be reinterpreted in 2D or 3D (or higher di- 
mensions) as a limit of a diffusion process where the bound- 
ary or initial conditions are replaced by fluid catalysts. It 
has been shown that pre-computing the diffusion process 
for an elementary catalyst case as a fundamental block of 
a class of differential equations, the computation efficiency 
can be greatly improved. Here, we explain how the diffusion 
point of view can be applied to decompose the fluid diffu- 
sion process per direction and how to handle non-symmetric 
discrete equations. The method can be applied on the class 
of problems that can be addressed by the Gauss-Seidel iter- 
ation, based on the linear approximation of the differential 
equations. 
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1. INTRODUCTION 

The iterative methods to solve differential equations based 
on the linear approximation are very well studied approaches 
[13], [2], [15], [3], [IT], 0], [16]. The approach we propose 
here (D-iteration) is a new approach initially applied to nu- 
merically solve the eigenvector of the PageRank type equa- 
tion [10], [9], [8], [6], 0, [UJ. 

The D-iteration, as diffusion based iteration, is an iter- 
ation method that can be understood as a column-vector 
based iteration as opposed to a row-vector based approach. 
Jacobi and Gauss-Seidel iterations are good examples of 
row-vector based iteration schemes. While our approach 
can be associated to the diffusion vision, the existing ones 
can be associated to the collection vision. 

In this paper, we are interested in the numerical solution 
for linear equation: 



X 



P.X + B 



(1) 



where P and B are the matrix and vector associated to 
the linear approximation of differential equations with initial 
conditions or boundary conditions. 

In [5], it has been shown how simple adaptations can 
make the diffusion approach an interesting candidate as an 
alternative iterative scheme to numerically solve differen- 
tial equations. In [12], a new approach based on the pre- 
computation of the elementary diffusion limit has been pro- 
posed. In this paper, we study the case of non-symmetric 
(from the diffusion point of view) linear equation and how 
we can very simply decompose the iteration method per di- 
rection for an improved convergence speed. 



2. GENERAL EQUATION IN 2D 

We consider the general linear (afhne) equation of the 
form: 



U(n, m) = a(+l, 0)U(n - 1, m) + a(0, +l)U(n, m - 1) 

(2) 

+ a(-l,0)U(n+ l,m) + a(0, -l)U(n, m + 1) + f(n,m) 

(3) 

for (n, m) £ Q, and with boundary condition: H (n, m) — 
g(n,m) for (n,m) £ dQ. We require that dQ, includes at 
least all the boundary positions (frontier) of fi. 

Let's call this problem DD2D (discrete differential 2D) 
equation problem. We define: \a\ = \a(— 1, 0)| + \a(0, — 1)| + 
K+l,0)| + |a(0,+l)|. 

Note that the approach proposed here can be directly ex- 
tended to higher dimension and also for neighbour positions 
in Equation ([2]) ((n — 1, m) etc) that are more general (it is 
just required that they are regular on the domain where we 
iterate the equations). 



Theorem 1 (Stability condition). If \a\ < 1, then 
DD2D has a unique solution. The solution can be obtained 
from the iteration of Equation ([2]) . 

Proof. The proof is straightforward noticing that the 
matrix associated to DD2D has a spectral radius strictly 
less than 1. 

Now, we associate to DD2D, the diffusion approach: the 
associated diffusion approach consists in (iteratively) apply- 



ing the elementary diffusion operation on (n, m) by: 
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F(n,m + 1] 


+ = a(0, +l)F(n 
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F(n — 1, m 


+ = a(-l,0)F(n 
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F(n, m — 1] 


+ = a(0,-l)F(n 


m), 
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F(n, m 


= 0. 




(9) 
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Figure 1: Collection/Diffusion approach. 

We first build the initial condition for the diffusion pro- 
cess by applying the elementary diffusion operation to all 
boundary position in dQ,. 

The DD2D can be equivalently written under the form: 

X = P.X + B 

where the initial fluid B on all positions of is defined by 
the superposition (sum) of / and those coming from the 
boundary positions. 

We recall that the diffusion approach (D-iteration) con- 
sists in iteratively solve the above equation with the fluid 
vector F and the history vector H (cf. [9]), using Q. 

Theorem 2. Under the stability condition, the limit of 
the diffusion scheme defined above converges to the solution 
of DD2D. 

Proof, ft has been shown in [11] that our history vector 
H corresponds exactly to U in |(2J when starting from H = 
(0, ..,0) and when applied for the same sequence of vector 
entries (than the diffusion process on the fluid vector F Q, 
one vector entry corresponding here to a spatial position 
(n,m)). 

Below, we show through simple examples how our ap- 
proach works very simply and how it can help us decompos- 
ing the diffusion process for a better efficiency. 

2.1 Catalyst limit in ID 

The general case of the linear DD1D operator associated 
to the diffusion is: 



H(n)+ = F(n), 
F(n + 1)+ = a(+)F(n), 
F(n-l)+ = a(-)F(n), 

F(n) = 0. 



(10) 

(11) 
(12) 
(13) 



Its elementary catalyst limit (cf. [12]) (j> is associated to the 
solution of (|10p with the initial condition g(0) = 1 and with 
the constraint that the diffusion is applied once at position 
and then the position behaves as a black hole (diffusion 
only on n / and with at +oo). The solution is here 



simple and explicit (for instance for a(+) > and a(— ) > 
and | a | < 1): we can solve the equation 

a(-)x 2 - x + a(+) = 

for the solution on 7L + , which is (f>°^(n) = r. 

1- v /l-4a(+)a(-) 
2aR 

and for 7Z~ , 0~(-n) = rl with 



with 



1- v /l-4a(+)a(-) 
2a(+) 

One can easily check that the above solution is the right 
one for all cases of |a| < 1. 

If we put a finite boundary condition g(N) = 0, we find 
the exact limit by successive compensations of the surplus 
or deficit of fluid at the two boundary positions: 
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. \—n) = r 



l-(r+r-) N ' 

l-(r+r-) N - n 



l-(r+r-) N ' 
In particular, if r = r + = r_, we have on [-N, N]: 

1 _ r 3(JV— 



<p (n) = r' 
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Note that <f> is the elementary function we can use when 
injecting fluid from the boundary position. When the value 
at the origin is not imposed (0 is not a black hole), we have 
to use the normalized version of cp; 



V ' l-a(-Ml)-a(+M-l) 
on TL without bound and on [-N, N], 

4> N {n) 



<j> N (n) 



l-a(-)ct> N (l)-a(+)(t> N (-l)' 



Remark 1. 0(0) = 1 by definition. 0(0) > 1 repre- 
sents the total fluid that comes back to 0. a(—)(j)(l) repre- 
sents the fluid that goes to the black hole coming from Zj + . 
a(+)4>(— 1) represents the fluid that goes to the black hole 
coming from 2Z~ . 

2.2 Example of ID differential equation with 
time dimension 

Let's consider a concrete case (cf. Q~]) of heat equation 
evolution in time in ID: 

d t U(t,x) =d x 2 U(t,x),(t,x) € [0,T] x [0,1] 

with initial condition U(0,x) = Uo(x) = sin(nx) (pre- 
heated metal stick) and boundary condition U (t, 0) = U(t, 1) 
(imposed temperature at the extreme points of the metal). 
Then we can discretize the above equation as (usually called 
implicit equation): 

U(t, n) - U{t - 1, n) _ U(t, n + 1) + U(t, n - 1) - 2U(t, n) 



At 



Ax 2 



(14) 



which can be written as: 



u(t > n) = TT2k U{t - 1 ' n ^ + TT2k {u{t > 71 + 1} + u{t ' n ~ 

First it is well known that the above scheme is always 
stable (cf. Q3]). Here, we have a(-l,0) = 0, a(0, +1) = 
a(0, —1) = f^ii an d Q (+l> 0) = 1 + 2 k • ^ ne initial condition 
is build by injecting Uo(x) to F(l,x): 

F(l,x)--a(+l,0)Uo(x). (15) 

Now, thanks to the freedom of the order in which we apply 
the diffusion (on the position choice), we do the following: 

• we first advance on the time axis once using only the 
diffusion with a(+l,0) (this is the application of (|15[) 
for time t = 1); 

• then we freeze the diffusion on time axis and diffuse 
only on x-axis; since we know the exact limit of the 
elementary catalyst solution (diffusion of 1 from x = 
with boundary condition equal to at N, which is 
4> N (n) (cf. SectionEjJ with r = (1 - y/l -4a 2 )/(2a) 
(a = a(0, +1)) we can diffuse all fluid for x — 1 to x — 
L x — 1 using the pre-computed elementary solution, 
and obtain the exact directional (x-axis) diffusion limit 
(without the need to compensate the surplus fluid at 
the boundary x — 0, x — L x ); note that we could 
also use ^(n) = r n /(l - 2rk/(l + 2k)), then decide 
to compensate the boundary values using iteratively 
<^°°(n) = r n or using once <j) N (n)); 

• then we restart the process for the next time t + 1 by 
diffusing H(t, ) with a(+l,0): 

F(t + l,x) := a(+l,0)H(t,x). 

Note that there is no approximation in the above scheme. 
We show the comparison of our approach (Dl) to the naive 
iteration of Equation (|14[) (GS as Gauss-Seidel) . 



Dl: 10s 




Figure 2: D-iteration: 10s (exact limit). 

2.3 Example of ID differential equation of or- 
der 2 

Let's consider the general second order linear differential 
equation: 

y"(x) + ay'(x) + Py(x) = f(x). 
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Figure 3: Gauss-Seidel: 37m. 




Figure 4: Gauss-Seidel: 2h30m. Not yet converged. 

Its discretized equation is: 

yjn + 1) + yjn - 1) - 2y(n) ^ y(n) - y(n - 1) 

Ax 2 Ax 1 ' 

+ Py(n) = f(n). (17) 

which can be written as: 

y(n) =a(-)y(n + 1) + a{+)y{n - 1) - 7 /(n), (18) 

with q(-) = ^— L—^ | Q ( + ) = _=gg« s _ a and 7 = 

Ai 2 

A sufficient stability condition for this equation is: A < 
min(^i, l/\a\, We can then apply the method as 

in the previous section using the elementary catalyst limit 
associated to a(+) and a(— ) (cf. Section [2TT} ■ 

3. CONCLUSION 

In this paper we showed that thanks to the diffusion point 
of view we can efficiently solve the linear equations asso- 
ciated to non-symmetric operators and that we could also 
exploit the idea of diffusion per direction for a faster compu- 
tation. This last idea of diffusion per direction is a promis- 
ing approach in the context of linear operators in higher 
dimension when the naive iteration method becomes really 
painful. Further exploitation of this will be addressed in a 
future paper. 
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